%% The script uses Jicamarca ISR data and ACE IEF data and tries to compute the
% transfer function between them.
% The jicamarca isr data is processed by the script file
% jcamarca_isr_processing.m
% results. I find that the tf is too noisy
clear;
load  /data/backup/mnair/ace/jcamarca_isr_fejer.mat eef min_time_length lt_start lt_end data_index;

load /data/backup/mnair/longp/OMNI_ELEC_new;

%load c:\manoj\projects\ace\Em.mat; %merging e field = -Vsw * Bt * sin^2(theta/2))
load /data/backup/mnair/longp/aplist.mat; 
fday_ap = fday_ap - datenum(2000,1,1);

% Initializations
%w=w_climate; %This copies the Julia drift (measured)-Julia drift (modeled) to w
ace_fday = floor(ace_all(:,1));
ap_lower_limit = 0; %lower limt if Ap
N_seg = 1;%the increasing counter for array JULI_SEG & ACE_SEG
N_data = 0;
phase_delay = 17;%minutes
mjd_date = datenum(2000,1,1);
plc = 'b';%plot color
tol = 1.8/24; % max time interval tolerated
%min_time_length = 32/24; % minimum time length required in decimal days
np = 0;
%lt_start = 11; % Ideal 9-15 (data permits 7-17)
%lt_end = 13;
nd = 1;
plc = 'r';
coh_plot_lim = 0.01;
ace_interp_method = 'spline';

des_int = 0.5;%desired sampling interval in hours
len = floor( min_time_length *24 / des_int ); % This number is the time length (also data length) in hours
% For changing this, first run the data selection in
% chamP-eej_spectra with desired min_time_length and save the
% results to  /data/backup/mnair/longp/eef_data_mod.mat

ace_gap_tol = ( min_time_length + 6/24 ) * 0.1 ;
ace_gap_len_tol = floor (ace_gap_tol*1440/5 );
% the gap tolerence is 10% of the ace data length

%mjd_limit =
% processing of ACE data
%down sampling and resampling to original interval

%The following script is to resample the ace data at a lower rate.
% The matlab resample apply a low-pass filter before resampling
% this allows for some remedy to the aliasing issues. This is a very
% important step. With out this sever aliasing issues can crrep into the
% spectra analysis.

L = isnan(ace_all(:,2));
y = interp1(ace_all(~L,1),ace_all(~L,2),ace_all(:,1),ace_interp_method);
%ace_all(L,2) = 0;
%ace_down = resample(y,1,des_int*60/5);
%ace_inter = interp1(ace_all(1:des_int*60/5:end,1),ace_down,ace_all(:,1));
ace_all(:,2) = y;
%ace_all(L,2) = 0;%making sure that non data parts are zero

%ace_all= [ace_all(1:des_int*60/5:end,1) ace_down];

% Arrange the CHAMP data into segments
% the climatology is limted to > 7 and < 17 LT

% iterations
CHAMP_SEG = [];
ACE_SEG = [];
TIME_SEG = [];

for i = 1: size(data_index,1),
    
    
    L = ace_all(:,1) + phase_delay/(60*24) >= eef(data_index(i,1),1) - 3/24 ...
    & ace_all(:,1) + phase_delay/(60*24) <= eef(data_index(i,2),1)+ 3/24 ;
    
    if sum(L) > 0,
        ace_time = ace_all(L,1) +  phase_delay/(60*24) ;
        ace_data = ace_all(L,2);
        % imf_bz_data = ace_gse_bz(L);
        %mean_imf_bz = nanmean(ace_gse_bz(L));
        %  mean_sw= nanmean(sw(L,1));
        L = isnan(ace_data);
        if any(ace_data),
            if  max(diff([ace_time(1) ; ace_time(~L) ; ace_time(end)] )) < ace_gap_tol && ...
                    sum(L) < ace_gap_len_tol,
                
%                 if ( sum(L) ~= 0 ),
%                     
%                     if L(1) ,
%                         X = find(L==0);
%                         ace_data(1) = ace_data(X(1));
%                     end;
%                     if L(end),
%                         X = find(L==0);
%                         ace_data(end) = ace_data(X(end));
%                     end;
%                     
%                     ace_data = interp1(ace_time(~L),ace_data(~L),ace_time);
%                 end;
%                 
                
%                 ace_down  = resample(ace_data,1,des_int*60/5);
%                 ace_data = interp1(ace_time(1:des_int*60/5:end,1),ace_down,ace_time);
%                 
%                 L = isnan(ace_data);
%                 ace_data(L) = [];
%                 ace_time(L) = [];
                %
                
                %  fprintf('Number of missing points on fday %d = %d\n', Julia_W(i).fday, sum(L));
                %if abs(mean(diff(ace_time))-0.0028) <= 2.7778e-005, % Use this with Burke data
                if abs(mean(diff(ace_time))-0.0035) <= 1e-004, %Use this with ACE 5 min averages
                    %  if abs(mean(diff(ace_time))-2/24) <= 1e-004, %Use this with ACE 2 hours averages
                    dummy = eef(data_index(i,1):data_index(i,2), 6 ) - eef(data_index(i,1):data_index(i,2), 15 )/1e3;
                    time_axis_desired = eef(data_index(i,1), 1 ):des_int/24:eef(data_index(i,2), 1 );
                    CHAMP_EEF = interp1(eef(data_index(i,1):data_index(i,2),1), dummy,time_axis_desired);
                    
                    ACED = interp1(ace_time, ace_data, time_axis_desired);
                    
                    L = isnan(ACED);
                    ACED(L) = [];
                    
                    L = fday_ap >= eef(data_index(i,1), 1 )& fday_ap <= eef(data_index(i,2), 1 );
                    mean_ap = mean(ap(L));
                    
                    if sum(isnan(ACED)) < 1 & mean_ap >= ap_lower_limit & length(ACED) >= len & length(CHAMP_EEF) >= len & mean_ap <=20,
                        
                        
                        CHAMP_SEG(N_seg:N_seg+(len-1)) = CHAMP_EEF(1:len);
                        ACE_SEG(N_seg:N_seg+(len-1)) = ACED(1:len);
                        TIME_SEG(N_seg:N_seg+(len-1)) = time_axis_desired(1:len);
                        % IMF_BZ_SEG(N_seg:N_seg+71) = IMFBZ(1:72);
                        N_seg = N_seg+len;
                        N_data = N_data+1;
                        
                        
                    end;
                else,
                    %fprintf('Day %d has some missing time stamp\n', Julia_W(i).fday);
                end;
            end;
        end;
    end;
end;

fprintf('Done ! N_data = %3d , Length = %d hours, LT = %d-%d, AP limit = %d color = %s, interp = %s\n', ...
    N_data, min_time_length*24, lt_start,lt_end, ap_lower_limit, plc, ace_interp_method);

% Coherence, phase and tranfer function
%JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
figure(1);
[Cxy_long,F] = mscohere(ACE_SEG, CHAMP_SEG,hanning(len),0,len,1/(des_int*3600)); %1/(5*60) = sampling frequency in Hz )
[Pxy,F] = cpsd(CHAMP_SEG,ACE_SEG,hanning(len),0,len,1/(des_int*3600)); %
[Pxx,F] = pwelch(CHAMP_SEG,hanning(len),0,len,1/(des_int*3600));
[Pyy,F] = pwelch(ACE_SEG,hanning(len),0,len,1/(des_int*3600));
phase = angle(Pxy);
[Txy_long,F_long] = tfestimate(ACE_SEG,CHAMP_SEG, hanning(len),0,len,1/(des_int*3600));
%tf = conj(Txy');
L = Cxy_long > coh_plot_lim;

Err_long = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_long)./Cxy_long ) ) .* abs(Txy_long);

Err_long_lg = 0.434 * Err_long ./ abs (Txy_long);
%errorbar(log10((1./F_long(2:end))./(3600)), log10(abs(Txy_long(2:end))),log10(Err_long(2:end))/2,log10(Err_long(2:end))/2,'g');
errorbar(log10((1./F_long(L))./(3600)), log10(abs(Txy_long(L))),(Err_long_lg(L)),(Err_long_lg(L)),plc, 'linewidth',2);
hold on;
axis([ -inf   inf   -3   -1]);

for i = 2:length(Cxy_long),
    fprintf('%06.3f ',(1./F_long(i))./(3600));
end;
fprintf('\n');
for i = 2:length(Cxy_long),
    fprintf('%06.3f ', Cxy_long(i));
end;
fprintf('\n');
%% JULIA data 
load /data/backup/mnair/longp/alldays JULI_SEG ACE_SEG;
JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
[Txy_short,F_short] = tfestimate(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60)); %Txy is agains calculated just to get
[Cxy_short,F_short] = mscohere(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60));
N_data = length(ACE_SEG) / 72;
Err_short = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_short)./Cxy_short ) ) .* abs(Txy_short);
Err_short_lg = 0.434 * Err_short ./ abs (Txy_short);
%errorbar(log10((1./F_short(2:end))./(3600)), log10(abs(Txy_short(2:end))),log10(Err_short(2:end))/2,log10(Err_short(2:end))/2,'g');

%errorbar(log10((1./F_short(2:end))./(3600)), log10(abs(Txy_short(2:end))),Err_short_lg(2:end)/2,Err_short_lg(2:end)/2,'k');
errorbar(log10((1./F_short(2:end))./(3600)), log10(abs(Txy_short(2:end))),(Err_short_lg(2:end)),(Err_short_lg(2:end)),'b');


